/*
 * DOS.h
 *
 *  Created on: 2011-1-8
 *      Author: zhang
 */

#include <math.h>

//create dos graph data
//template<class T>
//void createDoS(T* h_U,int M,long hits,char* fileName);

#define KPI 3.14159

namespace DOS{

	template<class T>
	inline T* createKernel(unsigned int m)
	{
		T* mu = (T*)malloc(sizeof(T) * m);
		for(int j=0;j<m;j++){
			mu[j] = (T)1/(m + 1)*((m-j+1)*cos(KPI*j/(m+1))
										   + sin(KPI*j/(m+1))*cos(KPI/(m+1))/sin(KPI/(m+1)));
		}
		return mu;
	}

	/******************************
	 * this function create a the dos just for a single point.
	 * mu: chebyshev moments array
	 * m: number of chebyshev moments
	 * x: the point at which the DOS value would be calculated.
	 ****************************/
	template<class T>
	inline T createDOSSingle(T* mu, unsigned int m, T* kernel, T x)
	{
		unsigned int j=0;
	  	T Re=0, T0, T1, T2;

	    //j = 0;
	    T0 = 1;
	    Re += mu[0] * kernel[0];

	    //j = 1;
	    T1 = x;
	 	Re += 2 * mu[1] * T1 * kernel[1];

	    for(j=2;j<m;j++)
	    {
	      T2 = 2 * x * T1 - T0;

	      Re += 2 * mu[j] * T2 * kernel[j];

	      T0=T1;
	      T1=T2;
	    }
	    Re *= 1/(3.14*sqrt(1-x*x));
	    return Re;
	}

	/*****************************
	 *
	 * this function create dos function
	 * mu: chebyshev moments array
	 * m: number of chebyshev moments
	 * s: number of samples selected for calculating DOS
	 *
	 *****************************/
	template<class T>
	T* createDOS(T* mu, unsigned int m, int s)
	{
		  T* kernel = DOS::createKernel<T>(m);

		  T* dos = (T*)malloc(sizeof(T) * s);
		  memset(dos,0,sizeof(T) * s);

		  int j;

		  for(int e=0;e<s;e+=1)
		  {
			T x=(T) (e-s/2) / (s+1) * 2;
			dos[e]=DOS::createDOSSingle<T>(mu, m, kernel, x);
		  }
		  return dos;
	}
};








